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This progress report describes work accomplished during 
the first six months of the current contract period (May 23 - 
November 22 , 1979). Progress has boon made on the following 
topics, for each of which a technical report is appended to 
this report: 

B.l Current Progress with the Numerical Model of Planet 
Growth 

B.2 Experimental Studies of Collisions 

B.3 Hydrodynamic Code Collision Calculations 

B.4 Orbital Element Changes Due to Gravitational 
Scattering 

B.5 Validauion of Mass-Shifting Algorithm 

B.6 Inclination-Type Resonance in a Dissipative Medium 

B.7 Theory of Rotations 

B.8 Origin of Asteroids: Velocity Pumping by Large 

Jupiter-Scattered Planetesimals 

Most of the above work was accomplished during the period May - 
July. 

More recently, Hartmann has developed a methodology to 
evaluate the rate at which megai-egolith increases its depth d as 
a function of total accumulate number of impacts on an initially 
smooth, coherent surface. This methodology involves tabulating 
the percentage of surface covered by craters penetrating deeper 
than depth d. This relation was then combined with the time 
behavior of early lunar cratering (Hartmann, ^ / Basaltic 

Volcanism Study Project, Chapter 8, in preparation) to give the 
megaregoli+-h depth on any rock surface formed (i.e., solidified) 
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at time t during the early history of the Moon. This rate 
of destruction of rocky crust can bo compared with rates of 
formation of rocky crust by solidification of the magma 
ocean, to show how early crustal structure evolved on plant 
surface ' as a result of accretionary cratering (see Figure 1) . 

Results of Hartmann's work were presented at the Lunar 
Highlands Conference, November 1979 (See Appendix Q.l) and a 
paper is being prepared for those conference proceedings. 
Recent results on the planet building process were also pre- 
sented (by Greenberg) at the DPS meeting in St. Louis (see 
Appendix Q . 2) . 


FIGURE 1: Three staqes of lunar surface evolution identified 

from quantitative models of megolith formation rate and magma 

ocean pooling rate. In the chaotic pre-crustal stage, the 

cratering is intense enough to penetrate ttrough lithosphere 

and mix ejected fragments v/ith molten magma. In the second, 

crust-forming stage, upper layers of crust are pulverized into 

megaregolith (resetting isotopic clocks?) , but igneous cumulate 

crust can begin to form at a depth sufficient to be shielded 

from cratering. This is likely to be tens of kilometers deep. 

Flows on surface daring this stage (4.3 to 4.0 Gy ago?) would 

be pulverized. In final stage, cratering has dropped to rate 

2 

where surface flows of plausible thickness (few 10 m) would 
not be completely destroyed. Their lower layers would thus 
serve as a source for rock samples. This stage (hence datable 
rock samples) would be younger than 4.0 Gy old. 
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RECENT RESULTS AND CURRENT RESEARCH 


Con! ents 

L.l Current Progress with the Numerical Model of Planet 
Growth 

B.2 Experimental Studies of Collisions 

B.3 Hydrodynamic Code Collision Calculations 

B.4 Orbital Element Changes Due to Gravitational Scattering 

B.5 Validation of Mass-Shifting Algorithm 

B.6 Inclination-Type Resonance in a Dissipative Medium 

B.7 Theory of Rotations 

Origin of Asteroids: Velocity Pumping by Large 

Jupiter-Scattered Planetesimals 
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APPENDIX B.l 


CURRENT PROGRESS WITH THE NUMERICAL MODE 


OF PLANET GROIVTH 


B.1.1 


APPENDIX B.l 


CURRENT PROGRESS WITH THE NUMERICAL MODEL OF PLANET GROWTH 


Our numerical model serves as a central focus for our 
multi-disciplinary study of planet growth. As we make 
progress in the various technical isr-ues related to the main 
problem, we continue to use those results to update and 
improve our numerical simulation. The model provides 
insight into the interactions of various processes which 
are likely to have operated at that time. 

The major new mechanism that has been incorporated 
into our simulation in the past few months is gas drag. 
Before our original simulation was constructed, one might 
have reasonably imagined that most primaeval gas had been 
removed before the epoch of collisional accretion began. 
However, our results (cf. Greenberg et aJ^. , Icarus 35 , 1, 
1978) have shown that substantial growth may have occurred 
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in only a few 10 yr. It thus became less easy to believe 
that gas had been removed before this growth, and it became 
important that we incorporate gas drag. From the analysis 
of Adachi et al. ( Prog. Theor. Phys . 56 , 1756, 1976) and 


following the notation of Weidenschilling (MNRAS 180 , 57, 1977) : 


B.1.2 


de ^ (2q/; g) e 
dt t 

e 


.770 + .6-U 


2g 


li = (g/’g ) 

dt t 


1 I .77e + .85i - 


2g 


da 

dt 


■2a 


97e -1- .64i - 


2g 


where 

■ nRT ^ - 

Ag = - < 0 

^ yr 

r = distance from sun 
P a r ^ with n a- 3 . 5 
y E molecular wt. of gas ^ 2.25 
R E gas constant 
T E temp. 'V 300°K 

g = solar gravitational acceleration 

also t^ is the e-folding time for velocity change due to 

gas drag. For large Reynolds numbers {applicable if planetesimals 

are ^ 10 m in diameter) , 

4- = 3Pp 

e PgAV 
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where 


D 

P 

AV 


planetcsimal diameter 

planetesimal density 

“9 3 

gas density 10 gm/cm in our nominal model 

Y _ Y = - ( — ^ V 

Kepler gas \ 2g / Kepler. 


Thus the gas acts to damp random motion (c and i). But, 
by .introducing a radial drift (da/dt) that varies with 
body size, it introduces a new source of relative '^'elocity. 
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Values of da/dt for our nominal gas density range from 10 
cm/soc for a 10 m :, jdy to only 10 cm/sec and less for a km 
or larger body. The numerical simulation program has been 
modified to include this differential radial drift in the 
computation of relative velocities. The relative velocities 
in turn may govern both impact velocities and gravitational 
stirring. In principal, gravitational interactions can 
convert radial drift into increased random motion. However, 
so far in application to our program (cases to be discussed 
in detail below) , conditions have been such that gas drag has 
primarily damped random motion, rather than stirred it. 

A great deal of effort has been expended in making a 
number of other modifications to our code that improve our 
confidence in its validity. We find that such changes have 
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had noglicjible cffoctr, on .previously published results. 
Nevertheles'^ , continued checking and validation of our code 
is important ’■ f we are to ensure that it in fact simulates 
the physical processes that it is supposed to model. 

One such change was required when we found irregularities 
in the behavior of the inclination of the largest body. The 
change in each time step, though small, was the difference 
between a large gravitational stirring increase and a large 
collision damping effect. For this and other reasons, a 
numerical procedure had to bo developed to avoid loss of 
precision due to round-off errors. 

Another problem was that under certain conditions, 
we found unaccounted-for changes in total system mass of 
a few percent. After considerable effort, we found and 
corrected several minor "bookkeeping" errors, so now mass is 
conserved to 8 figures. 

We also modified our algorithm for shifting bodies into 
an adjacent size bin as their masses gradually change due to 
cratering erosion or accretion. Previously, we computed the 
number of particles shifted out of the bin and added an 
appropriate number to the adjacent bin so as to conserve mass. 
We now compute the total mass of particles shifted and adjust 


numbers in both bins to account for the mass shift. This is 
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more meaningful physically, although both methods would bo 
equivalent for infinitesimal bin-widths. Use of the new 
algorithm does not change any previously reported results, 
but it gives us greater confidence in the continuation of 
growth to larger bodies to be discussed below. See also 
Appendix B.5 for discussion of continuing work in the 
validation of our numerical methods. 

As we shall show, the general conclusions of our earlier 

work (cf. Greenberg ^ / Icarus 35 , 1, 1978) still hold 

given the program modifications discussed above. Briefly, 

those conclusions were that 500 km bodies can grow from a 

. 4 

swarm of km bodies in only '\/2 x 10 yr. At the time they 
grow, most of the mass remains in km size bodies. In 
accordance with the theory of Safronov (1972, NASA TT F-677) , 
the random relative velocities remain small, on the order 
of km bodies' escape velocities. We stopped our runs at this 
point because we felt random eccentricities were so 
small compared to the "Keplerian" accretion annulus of the 
largest body that our computation of collision probabilities 
was invalid. Now we are extending our simulations forward in 
time to investigate evolution in a later stage of the growth 
process than had been considered before. A critical scientific 
question to be answered by extension to later stages (cf. 
Greenberg et , in Asteroids , T. Gehrels, ed., in press; 
Greenberg, Appendix A. 2) is whether an equilibrium condition 
of the type assumed by Safronov (most mass in largest bodies) 
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may have held during most of the growth period, or whether the 
"disequilibrium" situation, found by us in the early stage, 
persisted throughout much of the growth process. The answer 
to this question has implications regarding the size of the 
largest late-stage planetesimals and in turn regarding 
several observable properties of the current solar system 
(cf. Greenberg, Icarus , in press, and Appendix B.l). 

We now believe that we can extend our simulation past the 
epoch discussed in our earlier reports and that '.;e were too 
conservative in stopping our runs when we did. As Kaula 
(1979, preprint) has pointed out, although we found the 
largest bodies' orbits to be nearly circularized, the random 
motion of the smaller bodies, which dominated the population, 
was considerably greater. So our computation of the inter- 
action of the largest bodies with the bulk of the population 
was still valid. Furthermore, as we shall show, it remains 
valid because the velocities increase sufficiently as the 
larges bodies grow to several thousand km diameters. 

Now we examine the results of several recent numerical 
experiments. First consider an experiment with no gas, and with 
the "intermediate material" and initial conditions defined by 
Greenberg et (1978) (throughout this discussion the reader's 

familiarity with that paper will be assumed) . The size distri- 
bution evolution is shown in Figure 1. When the first 500 km 
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body is produood , occonl rio i t ios and inclinations of the 

-4 

bulk ot the population virt' a tew 10 , lame enough that we can 
continue the evolution as siunu'sted by Kaula's conuiient. As 
the ItU'iiest body iirows to almost 4000 km diameter, typical 
o's and i's are ^2 x 10 Wlii 1 e these increasetl V'l'locities 

accompany t lie grcnvth of the largi’st body, they are [irobably 
not: due to tiie I'ffi'cts of the lamest body. Ratlier, they 
are due to the shiftiiui of tlie mass peak in the sii:e 
distribution f rom 1 km towards larger values. Towards 
the end of the period stiulied, the largest body begins to 
contain vi substantial fraction of t lie total system mass, as 
it would in Bafronov''s equilibrium condition. However, it 
has already begun "runaway" growtli and is effectively detached 
from the continuous sice distribution, so it cannot effectively 
influence relative velocities (cf. Appendix A.l). Safronov's 
condition is not achieved at this point. Hy the time the 4000 
km body is produced, the smaller end of the distribution has be- 
come linear (a power law) . Moreover, relative velocities are so 
high tliat bodies up to i‘30 km are comminuting one another so 
that nearly half the mass of the system has been processed into 
small particles below the size range follow’ed in our numerical 
model. In reality, one might imagine a power law distribution 
extending down to such size (a crystal size or a few nun?) tliat 
particles' stremiths are larger than we have assumed. There 


B.1.9 


would bo a spiko in the distribution as the comminuted 
pieces pile up at that size. This material might experience 
substantial radial motion duo to drag (radiation or gas) 
and spiral out of the zone under study. But presumably 
it would be replaced by similar material moving in from an 
adjacent zone. (We will give more rigorous consideration 
to that idea.) The small particles would be an additional 
source of material for accretion by the largest bodies, so 
the runaway growth found in this experiment is probably an 
underestimate . 

There is a strong indication that Safronov's equilibrium 

condition may never bo reached. Implicit in his work, and 

explicit in recent publications (Safronov, 1979, preprint) 

is the idea that the largest bodies are all nearly the same 

size, as well as containing most of the mass of the system. 

Schematically, he assumes a size distribution (over the 

region which finally feeds a planet) to look like Figure 2 
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over a time scale 10 years: 



FIGURE 2 

In fact, we find a very different situation which 
can be summarized schematically as well: 



P3 anet 


n.] .11 


Most of tlio m.iss l^lay^3 ir* sniall btxiioa throuqhoul 
much of t he iirowtii, and may (,'Vcn ho concent ‘at od into dust 
particle sices. PI anetes imal s ramlc'i.n relative velocities 
stay very low compared to tlu' escape velocity of the larae.st. 
body. The larqcst bodies beqin to contain a substantial 
portion of the total mass only after runaway qrovs’th of the 
very larqest Iris bequn. There seems to be no prospect for 
Sa f ronov ' s equ i 1 ibr i urn cond i t i on to be reached . 

PI anetes imals destined to crow into full-siced 
planets sei'm to be selected by this runaway }>rocess within 

4 

a few 10 yr of our initial state. Conceivably, such 
multi-thousand-km bodies emerqed from tlie sice tli stribut ion , 
conserved their identities and evc>lved by a wide raneje of 
dynamical , qeophysical and qeochemical processes into tlie 
planets we know today. Alternat ivel y , a subsequent period 
of collisional and gravitational interactions amongst large 
bodies (cf. Cox et al. , Tejuapa 415, 1978) may have 

interrupted any continuity of identity from the first product 
of runaway through the final planets. A great deal more 
research will need to be done on late stage processes before 
such issues can be resolved. 

One potential shortcoming of our model is that the 
computed random velocities represent values averaged over all 
planetosimals of a given size whether or not they are in tlie 
"feeding zones" of the largest body. Presumably, in the 
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feeding zones velocities would bo considerably 
higher, perhaps giving local satisfaction of Safronov's 
assumption. Greenberg (1^79, Appendix A. 2) gives an 
argument against that possibility. Certainly, we cannot 
conclude that runaway growth always occurs, rather than 
Safronov's scenario, until the full range of possible para- 
meters and mechanisms has been explored. 

For example, consider how the numerical results are changed 

if we include gas as per the algorithm described above. The 

results of this experiment are shown in Figure 4. The 

basic results are quite similar to the previous case. 

Relative velocities of planetesimals are somev^7hat lower, 

yielding a 500 km body at t = 9000 yr. Characteristic 

-4 

e's and I's are about 10 at that time. By the time 
(5000 yi later) a 4000 km body is produced, velocities are 
an order of magnitude greater. Velocities remain low enough 
that only of the mass has been comminuted out of the 

system. Most of the mass is in the 4000 km body, and 
runaway growth seems to be underway. Velocities are low 
in accordance with Safronov's analytical work on relative 
velocities, for such a size distribution. But the low velocities 
and the runaway growth are in gross violation of the assumed 
conditions for his planet growth scenario. 

We have searched for a way to obtain Safronov's conditions. 
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For example, we note that Experiment No. 4 of Greenberg 
et al . (1978) , which uses parameters appropriate for a 
"weak material with low debris velocity," seems to be 
progressing toward a state consistent with Safronov's 
assumptions. At least the mass peak in the distribution is 
moving up in size almost as fast as the largest body 
grows, and velocities are considerably higher than for 
intermediate materials. We suspected that continuation 
of that experiment might lead to Safronov's condition. 

The results of continuing that experiment. Figure 5, while 
somewhat different in detail than the results of other 
experiments (Figs. 1 and 4) lead to the same general 
conclusions. The largest body reaches a size somewhat larger 
than 4000 km and e's and i's are ^^^.02. Despite the fact 
that relative velocities of planetesimals are high compared to 
velocities in earlier experiments, the comminution of significant 
mass out the small end of the distribution doesn't begin 
until after the largest body reaches 4000 km. This is because 
so much of the mass peak is shifted up to sizes which are 
immune to great comminution until after e and i become '^J.02. 

5 

After 4.5 X 10 yr about 1/3 of the mass has been comminuted 
out, but long before this, runaway growth has set in, with 
a "^6000 km (Mars-sized) body produced. 

We are also considering whether the addition of gas 
combined with varying material properties might give the 
Safronov equilibrium condition. .For example. Figure 6 
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shows a case with gas and the weak material of the previous 
case. Velocities are about half of what they were, the peak 
on the continuous portion of the distribution stays near 1 
km, and little comminution occurs. Runaway has 
begun. Figure 7 shows a case between that of Figures 5 and 
6: same material but less gas (t^ = 3 x previous case) . 

Results are similar. 

We conclude that, while we cannot now definitively 
answer the questions that we have posed about the later 
stages of planet growth, every indication from our work 
to date is that runav^?ay growth begins at an early stage and 
that Safronov's assumed equilibrium condition is never reached. 
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EXPERIMENTAL STUDIES OF COLLISIONS 
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EXPERIMENTAL STUDTBS OF COLLISIONS 
Amos Vert icq] Gu n 

Dr. Hartmann completed a highly successful 8-day 
run of experiments at the Ames Vertical Gun facility in 
Juno 1979. This was the first experimental run since the 
facility reopened after about two years of inactivity. 

Dr. Hartniann obtained collision and impact data in tlie 
hitherto-unstudiod velocity range around v = 100 to 150 m/s 
using an air gun designed by the Amos staff for those experi- 
ments. Data on collision outcomes are needed at these v's; 
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and on up to some. 10 m/s, where data were gathered by 
Gault et (1963) . 

Films and sample measurements were obtained at Ames 
which will give data on the following: 

• Ejecta velocities for basalt rock collisions at 
about 100 - 150 m/s. 

• Ejecta velocities for collisions of artificial 
rock aggregates at velocities about 10 to 150 m/s. 

• Ejecta velocities for material blasted out of 
powdery regoliths in vacuum at impact ■''clocities 
around 100 - 150 m/s 

• Size distributions of fragments in the first two cases. 
All of those data, when reduced, will be used for refining 

our computer modeling of collision events. The only publislied 
data on those subjects refers to collisions at less than 
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60 ni/s-., oi- Mi-iMti'i- tluin about 1 000 m,/u. Roduotion of l ho 
datii will 1 h' puriuioil tlurina L ho ourronl yoar. 

Throo oxamplor. of rouultr. will bo il i souiu-.od lu'ro to 
illuntrato prcHjro-'.M . Firrd, I’iuuro 1 i ;; a photoqrafdi 
showinq tho rosult of fi r i nq a .11 q proiootilo at 117 m/u 
into a bao.alt powdor (uomowhat ooarui'i' tlian lunar roqolitii; 
powdor :’>i;'o tl i ot r i but i on war. moarui'od by r. iovinq) in 
noar vacuum at mb prosruro. NoLo tho i-ay pattornr, rimilar 
in plan to tlioro on tho moon. Cartlboard rinqr woi'o urod to 
collect and woiqh ojocta at various annular di.-'tancor. I’ilmr 
of such I'Xporimont r. at about 7000 frames per socomi qivo 
moasuros of trajoctorios of o.iocta jMrticjos. Traioctory 
and mass moaiuiros will qivo volt^city d i s.t r i but i c'us of elect a. 
Those shoultl cTivi.'' tho i irst ri'liabli.' ilata on mass qain v«.'rsus. 
mass loss in modest -vol oc i ty collisions amonq p>l .uu't osinui 1 s. 
of different rises. Some pi'evious estimates of such phenomena 
relied on extrapol at ions of much hiqher-speed cc>llision d.Ua, 
such as pioneerinq data of (lault and others at around 6 km/sec. 

Second, Piqure 2 shows the si:u' di st ri but ions measured 
from two similar- collisions of 21 q and 166 q projectiles into 
rouqhly 3000 e artificial rock aqqreqate tarqets in which 
the typical mass of the constitiu'ut qrains was about 0.2 q. 

Tho figure shows that the mass distribution was about N " 

which is consistent with the distributions found durinq breakup 
of solid rock tarqets, showinq that the same f raqmentat i on 
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algorithms can bo used to represent these materials. The 
maximum in the size distribution falls at about the size 
of the original constituent grains, whose bonding into the 
aggregate had weaker strength than the grains themselves. 

This gives a handle on the mode of breakup of early planet- 
esimals that may be loose aggregates of small particles. 

Another feature of the graph is the relation of the original 
masses, largest fragments, and other fragments. This 
bears some superficial similarities to observed size distributions 
of small-body populations in the solar system today. 

Third, Figure 3 shows an updated approach to measuring 
and understanding strengths of planetesimal targets against 
disruption by impacts. Earlier experiments on which we based 
our estimates of "impact strengths" (e.g., Hartmann, Icarus 
33 , 50) were based on impacts into semi-infinite targets, in 
which only the projectile fragmented. The new experiments 
involved collisions with various mass ratios, where 
both target and projectile broke. These experiments tested 
the hypothesis applied in our numerical simulations: 

The energy is partitioned half into the target and half 
into the projectile, regardless of the target/projectile mass 
ratio. In Figure 3, the mass ratio of largest fragment to 
initial mass is plotted against the energy density delivered 
to the object under this hypothesis. It can be seen that 
for a wide range of conditions (projectile mass/target mass = 
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0 to 1; velocity = 10 to 150 m/s; data for projectile and 
data for target) , a rather well-defined curve divides the 
cratering regime from the regime of "catastrophic fragmentation" 
(the case in which the largest fragment £ half the mass of 
the original. The limiting energy density is rather well- 
defined near 3 x 10 ergs/g in these new data. This implies 
an "impact strength" of 'v^lO^ erg/cm^ . Since both the targets 
and projectiles follow the same curve, we conclude that our 
hypothesis of 1/2 of the energy being delivered to each is 
reasonable. This impact strength is about half of that 
found in previous interpretations of experimental data where 
all the energy was assumed delivered to the target. 

S.A.I. Santa Ana Facility 

A program of impact experiments was initiated during 
the past year in cooperation with the SAI impact laboratory in 
Santa Ana, California. The experiments performed to date were 
part of a no-cost precursor program to see if this facility 
was suitable for carrying out experiments relevant to 
problems of planet formation. The experiments carried out are 
summarized in Table 1. Targets were fine grained, irregular 
shaped basalts showing evidence of surface weathering 
collected on Sentinel Peak near Tucson. Sample No. 4 had 
several a priori hairline fractures, but the others appeared 


B.2.8 


K nj 

rH -H 

W (TJ +J 
W C -H 
(0 -H C 

s: M 


r^vxJcno OrotNoo 
•Tioocno criinfNCN 


C fd 

w s: 

I 

1 

1 VD 

VO 

in rr 

CO 


CO 

CO 

•p -p 

in 

m 

r»j 


CO 

in 

CO 

0 (U -N 

• 

• 

CN 

• • 

« 

• 

• 

• 

nj t7> CO 

iH 

r}< 

tN 

CN 

n 



CU Cn 
E fO D 
H E-( 



r- 


ro 


(N 

ro 

in 



• 

• 

• 


• 

• 

K 


o 



CTk 

o 

iH 

ro 

(N 

oc 

m 

in 

o 

O 

CN 

(T» 

oc 

00 

o 

m 

in 

o^ 

CT> 

00 


CO 


H 

m 










H 

■p 

in 

o 

o 

r- 

ro 

H 

cn 

o^ 


d) 

, 


t 

• 

• 

t 

• 

• 

< 

Cn to ^ 

LO 


00 

CT» 

cr^ 

CN 

CO 

in 

§ 

p to e 

CN 


H 

o 

o 

O' 

o 

cn 


(d (d CTi 


VD 

CTk 

0^ 


iH 

iH 

CN 


ei s 







n 


Eh 










;s 








o 


C 

■p ' 

<T> 

a\ 

CO 



00 

fH 

cn 

CO 

U to to 

O 


in 


o 


00 

in 


fd P to 

CN 


CN 

in 

CN 


o 

00 


a 0) P 



• 

■ 

• 

t 

« 

• 

o 

e c cu 

ID 

rH 



CN 


r-H 

■H 


H w 









>H 





















tu 









s 

H ^ 









s 

•IH p 









tJ 

■p & 



rH 

m 

o 

(N 

LO 

o 

CO 

0 


no 

fO 

in 

CO 


00 

CO 



tu 

CO 

in 

o 

CO 


CO 

m 




•n to 

• 

• 

• 

• 

• 

« 


• 

•P 

• • 

0 to 

CO 

CX3 

in 




in 


to 

H 

P Id 

H 








o 

W 

ft s 









ft 

PI 










•P 

CQ 

•p 









C/} 


u 'd — 









tu 

Eh 

(d (U to 

CO 

CO 

o 

CN 

cn 

o 

rH 

00 

tJl 


ft tu \ 

cn 

CN 

cn 

in 

o 


CN 

cn 

p 


E fte 

CN 

H 

r-f 


cn 

H 

CO 

CN 

(d 


fO 

^ (U 

4-1 m 

<D 

■pew 

CJ -H C 

rtJ -P (0 

&I (U 

e (u e 

■H O 

P C! 


^ q 
=«: O 


■f~■^—■f- 
iH -k -K * O 
nH*^DVo^oH'4*c» 


1 


B.2.9 


unfractured. 

Our results are plotted on Figure 3 for comparison with 
our Ames experiments. Note that the "impact strength" 
appears to be about a factor of 2 less than for the Ames 
experiments. This may be due to an intrinsic relative weakness 
in the field material. Based on the measured impact strength, 
these basalts would fall between the "intermediate material" 
and "solid rock" cases of Greenberg et aJ^. (1978 , Icarus 35 , 1) . 
Alternatively, it is conceivable that the glass projectiles 
deliver more than half the impact energy to the basalt targets. 
This will be tested by repeating similar experiments with 
basalt projectiles. Such results on the dependence of 
collision outcomes on material properties may have important 
implications for solar system evolution. 

Other aspects of the collisional model of Greenberg 
et al^. may be compared with these experimental results. For 
example, outcomes from our Santa Ana "cratering" experiments 
(i.e., those in which only a small portion of target 
material was removed) can be compared with the results pre- 
dicted by our numerical algorithm. For solid rock, the mass 

of ejecta is computed by "energy scaling" with a scaling 

-9 

coefficient, K = 10 gm/erg. However, we find that this 
predicts much smaller amounts of ejecta than the experiment 
yielded (Table 2) . 
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TABLE 2; PREDICTED VS. EXPERIMENTAL TOTAL EJECTA MASS 


TARGET 

ROCK 

PREDICTED EJECTA MASS 
(K = 10 ^ cgs) 

EXPERIMENTAL 

RESULT 

3 


2.6 gm 

95 gm 

6 

(first shot) 

0.2 gm 

8,3 gm 

6 

(second shot) 

0.03 gm 

0.4 gm 

6 

(final shot) 

1.1 gm 

99 gm 
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value K 4 X 10 (greater than that for craters formed 
in sand) would be required to describe such results. "You 
do experiments and you learn things," a famous scientist has 
noted. In this case, the discrepancy between our experimental 
results and our numerical algorithm (which was based on 
a compilation of previous experimental data for cratering) 
is due to a qualitative difference between the damage 
inflicted by these experiments and the formation of craters. 

In this case, the damage is of a chipping and/or spalling 
nature. Moreover, the irregular shapes of the targets have 
protrusions (relative to a spherical shape) which are more 
prone to breaking off. For these reasons, the simple "cratering" 
algorithm is just not applicable. 

We plan further experiments of this type in order to 
develop an algorithm to model these chipping events. At the 
same time, we need to consider both experimentally and by 
comparison with our knowledge of small-body properties, the 
conditions under which this algorithm should be applied. For 
example, asteroid observations show that irregular shapes are 
common over the km to '^^100 km size range. Thus, the planetesimals 
in our numerical simulation perhaps would be treated more 
realistically with a "chipping" algorithm under certain condi- 
tions. The continuing v;ork in this area allustrates how our 
numerical planet-growth simulation evolves in response to our 
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detailed technical studies of relevant processes. 

Next we consider the size distribution of debris created 
in the examples of debris of catastrophic fragmentation 
(Target rocl^s #4 and #8) . The size distribution predicted 
by the numerical algorithm are plotted for comparison with 
the experimental results. The impact strength inferred for 
these rocks from Figure 3 has been used to construct the 
predicted distribution. The predicted size of the largest 

fragment agrees well in each case with the experimental re- 
sults. However, the remainder of the fragmental distribution 

is somewhat different from that predicted by the Greenberg et al . 
algprithm; in particular, the predicted slope is steeper than 
is actually observed, as is shown in Figure 4. Further 
analysis is needed, both to extend the data base and to 
understand the cause of the discrepancy. 

The modeling of the large size end of the distribution 
may need to be re-examined. For example, Sample 8 produced 
three large fragments of approximately t.he same size rather 
than the predicted single object. The excess number of large 
objects mandates a lower slope in the 


• £xfCY'ih^0Hf [iesult dock 

Al 9or'ifkM Fr'fdictf^n Rock 
! _ 

X Expert htewf j^esul"^ Rock 
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remainder of the distribution relative to that theoretically 
predicted. The other sample fragmented into a distribution 
that is better represented by two power law segments, with 
the large size portion having a steeper slope than the small 
size end. In both cases, the smaller size population is 
significantly overestimated by our theoretical algorithm. 

SimiJar "non-linearities" in the size distribution of debris 
were found by Fujiwara ^ ^ (Icarus 31, 277, 1977) . 

The experimental program to date has provided useful 
results for validating our computational algorithms and pointing 
toward areas of deficiency in the models. The Santa Ana 
facility does appear suitable for carrying out further 
experiments there this year using their holographic imaging 
system to record the mass-velocity distribution of ejecta from 
high-speed impacts. 
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HYDRODYNAMIC CODE COLLISION CALCULATIONS 

Large computer codes developed to carry out hydrodynamic 
flow calculations have boon applied to study hypervelocity 
impact cratering by J.B. Dryan, D.E, Burton, M.E. Cunningham 
and L.A. Lettis, Jr. ( PLSC 9_, 1978, p. 3931-3964), J.D. O'Keefe 
and T. J. Ahrens ( Science 198 , 1249-1251, 1977; PLSC 7, 3007- 
3025, 1976), and R. Bjork 3379-3387, 1961), and other 

groups. Wo initiated a pilot project during the past year to 
try to extend hydrodynamic calculations to the case of 
catastrophic disruption of finite targets by high speed 
impacts. Our second objective was to obtain a capability of 
predicting outcomes for certain of our experimental shots and 
comparing the observed and predicted results. 

In cooperation with SAI's Huntsville office, a two- 
dimensional Eulerian finite difference code, HULL, was 
applied to invertigate tv^o problems. Computer time and labor 
to carry out those calculations wore contributed by SAI . The 
first case treated was a simulation of the Meteor Crater 
impact event for comparison with published solution of 
Bryan ejt al. (1978) , using the initial conditions of Bryan ejt al . 
The comparison is still underway (the donated effort is done 
on a time available basis) . However, readily available para- 
meters calculated by HULL are in good agreement with the pub- 
lished solution. For example. Figure 1 compared projectile 
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penetration depth as a function of time from HULL with the 
solution of Bryan o_t al . 

The second case investigated with HULL was a simulation 
of two 1 km radius granite spheres colliding at 1 km/sec. 
Figure 2 shov-/s a cross-section of the impacting bodies 
at several intervals based on HULL output. The 
calculation, such as velocity distribution of ejecta, will 
be carefully evaluated to see if they are plausible. 

Fragment sizes based on this calculation must bo treated 
with caution for material strength failure is the only 
fragmentation mechanism considered in the model . 
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ORBITAL ELEMENT CHANGES DUE TO GRAVITATIONAL SCATTERING 
Introduction 

Collisional velocities play a fundamental role in 
determining the evolution of a planetesimal swarm. If 
velocities become too large, collisional destruction and 
erosion results and the distribution proceeds to grind itself 
down, as is happening today in the asteroid belt. On the other 
hand, low encounter speeds result in accretion and growth of 
larger bodies. Collision speeds between particles are 
altered by interaction with other particles in the distri- 
bution - increased due to gravitati^.-nal scattering and 
decreased by collisional damping. The model for velocity 
changes due to gravitational scattering previously in 
the program was developed for initial stages of planetesimal 
growth as described by Greenberg ^ a_l. , (1978) . We are now 

extending and improving our model for velocity changes due to 

close encounters of bodies, in particular with application 

2 4 

to scattering by intermediate- and large-sized (10 - 10 km 

diameter) bodies. 

The basic problem is to find changes in mean orbit elements 
due to a large number of random gravitational encounters between 
n^^ bodies, each of mass m^, moving on orbits uniformly distri- 
buted within increments Aa, Ae and Ai centered on a^^, e^ and i^ 
and a second population r \2 described by another set of parameters. 
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The other orbit elements - node lonqitude, perizose longitude 
and mean anomaly are assumed to bo uniformly distributed 
thru 360°. The patched or sequential conic model is adopted 
to simulate what are really three body interactions. 

In this model, m^^ and m 2 are assumed to move on unperturbed 
conic orbits about the Sun when their mutual separation dis- 
tance exceeds a critical distance. Within the critical 
separation distance, they move on conic orbits about their 
mutual mass center. The critical distance is referred to 
on the sphere of influence (SOI) which depends on the masses 
of the bodies and their distance from the Sun. The most common 
criteria for the size of the SOI, the Tisserand criterion, 
leads to a SOI radius 



Comparison of results from patched conic models and 
fully integrated orbits (see Figs. 00 and 0) indicate 
generally rather good agreement except for distant, very low 
relative velocity encounters. By ignoring perturbations out- 
side the sphere of influence, the patched conic model under- 
estimates the frequency of small changes in orbital elements. 
The effect of a close gravitational interaction between two 
particles primarily produces a change in the heliocentric 
velocity vector of one or both particles. In the patched 
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conic model, the scattering process produces a rotation of 
the relative encounter velocity vector U but does not 
change its magnitude. It is this rotation that changes the 
heliocentric velocity vector. Provided the encounter time 
is small compared with the orbital period and the SOI is 
small compared with the orbit, the position vector is 
essentially unchanged during the encounter. Hence the post- 
encounter (scattered) velocity vector can be transformed 
to equivalent orbit elements a, e, and i. By averaging 
over all possible encounter orbits, the mean value of the 
scattered velocity vector may be found which can then be 
expressed in terms of mean scattered orbit elements. 

In the following sections we describe work in progress 
on two complementary methods of calculating changes in mean 
orbital elements due to close gravitational encounters based 
on the patched conic model. The first method (a) derives 
the probability distribution of orbit elements of the scattered 
body based on the mean deflection which is calculated with an 
approach analogous to that used by Wetherill (1967, JGR 72 , 

2429) in calculating orbital collision probabilities. The second 
second method (b) is to statistically average over all 
possible encounter orbits using the deflection angle probability 
distribution of Weidenschilling (1975, 145-153) and a 

uniform probability distribution of relative inclination angles. 
These two methods are used to determine the mean scattered 
heliocentric velocity vector from which the mean elements may 
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be found. Last year, reviewers of our proposal complained that 
we had not used enough equations to communicate our ideas. In 
this report a mathematical approach is appropriate and equations 
will be used liberally.' 

Method (a ) : 

One approach to the problem is to extend the method 
developed by Wetherill (1976, JGR ) for computing collision 
probabilities to the computation of probabilities of orbital 
element changes. First, we review Wetherill 's method 
using a slightly different, but equivalent, mathematical 
approach that will be readily extendable to the velocity dis- 
tribution problem. 

How much of the time are two different orbits of given 
a, e, and i within a distance t of one another? The "field 
particle" elements will always be denoted by subscript o. 

The orbital elements Q, to, F.to, and to^ (defined in Fig. 1) 
are assumed to vary uniformly with time so there is an equal 
probability of their having any value. 



Assume that the closest approach of the two orbits is in 
the vicinity of one of their mutual nodes. The position of a 
mutual node is given knowing Afi = and i and i^ . This 

geometry is shown in Figure 2 . 
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Note that Uq/ u and i* are woll-dofincd functions of i, 1 q 
and Afi: 

cos i' = cos i cos i = sin i sin i cos Afi (o<i'<180°) 

o o 


u = arc 
o 


tan sin An/(cot i sin i - cos A2 cor. i ) (2 values, 

*- o o^ 


180 


cos u = cos u cos AQ - sin u sin AH cos i 

O O Oi 


sin u = sin u sin i /sin i 
o o 


unique value, given u 
selection ^ 


Consider a coordinate system with origin where field orbit 
crosses line of mutual nodes, with X directed outward along 
line of nodes and XV plane in the plane of the field 
particle's orbit: 



Pt. B{X,Y,2) is on 
test orbit. Pt. E is 
its projection on Xl’ 
plane. 


apart) 
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Note <Xq is the angle the field orbit makes with the node 
line nndc<is the angle the test orbit makes with the node 
lino. For a Keplcrian orbit 


ct - 


e ■^| I - c^^(u-oS) 


I r e ' cv^ 

with a similar expression for oCq. Also 

^ 0 - g/) 

/ -f e c^^u.-cS) I i e^. 

The distance from any point B to the nearest point on the 
test orbit (assuming constant velocity, i.e., linear 
trajectories, near the mutual node) is d = whore 

z, ^ " i, oc>^ c ^ 

2 ^ oc /o-<^ I 

2 

Setting d (d )/dX=0 we find the value of X (i.e., the particular 
position on the test orbit) which minimizes ld|: 

(I 

6 r e <x /•Kn. i 

5 oi. ^e> ^ 


With this X, 


L - — 




, a 




2 — /(l-iTv o<.o 
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Hcncc, 






. ^ 




(i -x^y 


or, in a foi-m more symmetrical in terms of field and tost 
particle 

g t . 1. -t 

*J- 

^oC ^cCe> ^ ~ 2 

If dj^in is loss than t , the collision cross-section radius 
of the two bodies, the orbits pass close enough to one 
another that coliision is possible. Vv’hat can wo say about 
the dit;j^j^=0 surface in ( 0 , 00 ,^,^^) space? First note that 
the goomctx'y is independent of the absolute values of 
fi and 0.Q and only really depends on Afl. So wo want the 
dmin”0 surface in (oj,UQ,Afi) space. Setting aside the 
special case sin i'-O, dixiiirO only v/hore r.=0, i.e. where 

I + 

This defines a surface in (w/Ui^/Af!) space. 

Consider the cross-section of the dj^;j^n~0 surface at 
AS1=0. Note the linear relation between cos(u-t)) and cos(Uq-cOo). 
The ambiguity in arc cos and in u and u^ allows typically 
four values of 6J (if any) for which dmij^--0, given 
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all the other orbital elements (cf. Wothcrill's Fig. 1). 

For any other Ml, the locus has the same size and 

shape, except that its position is shifted according to 
the values of u and u^. We can consider an equivalent space 
(An, C, C^) where 



Fig. 4. 


An=constant plane 
a 


has d . -0 a 
min 

straight lino, d 
<T in a region 
defined by AC. 


min 


There is some region about the straight line 

which d . <T. How big is this region? Consider 
min 


Take 


^ cJ^hUh 

JcT" 




( + eC J 


1 f *] 

' O-cy 

L -yf <UXcr J 






Then 


AC -- r/(3cC:„/9c)l 

L-o 


r-* 
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AC = 

r 



[ V 



at 


! • ’ / 

w/iere 

Oi = 

e Vz-c’- 

u, h c*(. 




l-teC 


1 * e^Cc. 

at\cL C - 

[0 


aC/-e’) 

0.c(t-(^c) 



Therefore, AC is of form tf(C^). The differential area in 


which d . <T is 2tf (C ) dC . 
mxn o o 


The factor 2 accounts for 


both sides of the d . =0 line. We really want the area 

mxn 

in (u-£o, u^-(X)q) space, over which the a priori proba- 
bility distribution is uniform: 




An additional factor of 2 is required for the ambiguity 

of 2 mutual nodes, so the area for which d . <t is 

mxn 



■■ Such an integral, when extended over d(AJ^) will give 
the volume in (u,w^,AP) space in which • However, 

that is not a sufficient condition for collision. We also 
need to compute, for a given / 'the probability of 

collision, i.e., the probability that both bodies will be 
near where their orbits come together at the same time. 
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Consider the coordinate system shown in Fig. 5 in which 
the z axis is defined by the direction of d . from the 
tost to the field orbit and the y axis is parallel to the 
field body velocity U^. The velocities at closest approach 

-JV 

are assumed to be the same as at the mutual node. U lies 
in the xy plane. When the test body is at the 



origin, the field body is at y=); (t=o) . Separation of the 
two bodies is A: ^ 

Taking we find 

aL = y 

when 

t = ' U]ylu ^here Lid '' j U - , 
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For a collision wc require A . <t: 

^ nun 

, , = , 

roquiro )J < rZT 

{J C-T'O' "V^ 




If c/h\l^^ = whoro O'^K'^ljwe require 


Thus the probability of collision for any given is 

probability if Therefore the average probability 


IS 




times the probability if d . =^0. 

nun 

Now, given that what is the collision probability? 


Y 



Fig. 6: Same as Fig. 1 except 

d Field orbit is 

mxn 

in XY plane; Test orbit 
is out of plane. 


At t=0, test body is at point P with {:x.,y ,z)~ {,i\ cc^°^c> ^ C> ). 

Wetherill shows that collision occurs if 

I n I < r CLi / [ Du “ ( Udy y = 

Probability of collision each time test particle arrives 
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at oricjin is (timo during which / (orbi la] peritid) 


- £zk . J_ 
a X 


From Kcidcrian theory we know U : 

o 

1 1 ^ „ G H t e/ + 2 . (<■<■( ‘ ^0) 


and also 


T, = 




ZlT <xj 

'WW' * 


Probability of collision per time averacied over all d . <t: 

' nun 

JL l2i_ J. 

U, XT 

The not probability is that expression times the fraction of 
total (w, space for which x: 


P=i -3 




AJL->-0 


This probability has been evaluated for several cases for 
comparison with Wetherill's result. In particular, the 
hypothetical asteroid "Astrid" (a = 2.75, e = 0.2727, i = 0.2760) 
has been considered in pair with various real bodies. The 
following table compared the results we obtain with those 
obtained by Wetherill (method equivalent to ours) and earlier 
by Arnold (simplified method) . The probabilities are given 
for 
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TABLE: Comparison of three ways of computing probabilities 

of collision with "Astrid" 


Object 

Arnold 

Wetherill 


1948 EA 

1.92 

3.10 

2.49 

Apollo 

3.03 

4.22 

3.24 

Adonis 

4.95 

4.13 

3.92 

1950 DA 

2.19 

3.90 

3.13 

Geographos 

0 

0 

0 

Icarus 

0 

0 

0 

Encke 

2.56 

3.49 

2.91 

Brorsen 

0.79 

0.94 

0.81 

Grigg-Mellish 

0.020 

0.022 

0.021 


Given that the method described above was supposed to be 
equivalent to Wetherill's, why do the numerical results 
differ? We can suggest a couple of explanations, although 
there may be other, more important ones that we haven't thought 
of yet. First, Wetherill's assumed (second column, p. 2434) 
that the collision probability for a given 0^ H u^ - can 
be integrated uniformly over 0^ to get the average probability. 

In fact, some values of 0^ may be more probable than others 
because of the eccentric motion. Second Wetherill used an averaging 
system in which he evaluated the probability at a number 
of randomly selected values of 0^ and averaged them, rather 
than doing a numerical integration as we did. Whatever the 
reasons, our results agree to the level of precision required 
for Wetherill's 1967 conclusions. We should (and we intend 
to) be sure to understand the differences if we are to make 
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Other applications of our theory, however. 


In a gravitational interaction the relative velocity 
vector is rotated towards the direction of the minimum 
separation vector. Tlie previous discussion includes eval- 
uation of all the elements of the required geometry. The 
d . can be expressed in a vector form pointing from the 
field to the test particle's orbit (in XVZ coords, see Fig. 3) 



/"» . 


I L_ 

1 (^1 + ^ 1 

L, 



— T 


a 




This vector is the lino of minimum distance between the two 
orbits . 

The vector of minimum separation is given by considering 
the values of 7^ and of t at minimum separation as in the 
discussion of Fig. 5 and then rotating to the (>(,Y/Z) system 
of Fig. 3: 



.where t was found earlier and 0(in the expression for t) 
comes fromiX,^^^, and i' by spherical trig. 

Now the field particle's velocity relative to the test 
particle is y=U^-U. 




V 


> r 
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1 


^triin 


Fig. 7 


m 


This relative velocity is rotated by angle 2x where 

A 




^ v-'/i 


Q.^ (hn + 


vn- v^, 


a well-defined solution in tenns of given orbital elements. 
Let the components of Au in (x,y,z) be (6^ 6 „ S^) , These 

are the same components of AV used by Danby (1962) in his 
equations for the orbital element changes which correspond 
to a given AV. So, v;e get 




iTel 


VL(, O 






and similar changes for Ae^» etc. Thus for any given 

2 sets of orbital elements we have the change in elements 
at closest apprach near the nodes. Note that these equations 
work only if closest approach is near the nodes. Therefore 
we must set all changes equal to 0, if is greater than 

•some critical value such as the "sphere of influence". 


Also require that ea, ia, e^a , i a >> that critical distance, 
^ o o o o 
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Conclusion ; We can obtain the probability distribution for 
orbital element changes by computing all changes over uniformly 
space oj, AlT2 and mean anomaly of the field particle when 

the test particle is the point of By taking 

the average in this way, all changes will be "per orbital period 
of test particle". 

In calculating the changes in orbital elements, we 
really should consider both nodes. However, except in special 
cases, only 0 or 1 of them will have an interaction within the 
sphere of influence. 


Note ; One thing we should have shown earlier is how 

n is related to the mean anomaly Mq of field particle when 
test particle is at d . point: 



FIGURE 7: 

We know that the mean anomaly 
at the mutual node is given by the 
usual relation between mean anomaly 


and true anomaly, u^-u^: 

"’noae=* 


which you get from Kepler's equation. 

Then a time /^o when the field particle 

is at its d . point. At that time the field particle is 
min ^ 

a distance U (M --M , )/n past the node. Thus 

o o noae o ^ 
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rl= U Q 

where Q is the distance from the node to the d . point on 

min 

the field particle's orbit. 


What is Q? Refer to Fig. 3: 

Y 



where X and L have been 
evaluated previously. 
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Method (b) ; 

An alternative approach makes use of probability 
distributions for the rotation angle of the relative velocity 
vector and for the plane of the relative encounter trajectory. 

In the following analysis we use the equivalent single body 
scattering model. 

For simplicity, assume m2 >> m^^ so that the center of 
mass is at m2 (when the masses are comparable, the relative 
orbits can always be expressed as equivalent one body orbits) . 

At any position on the scattering body (m2) orbit where close 
approaches between m^^ and m2 are possible, we assume that 
there will be sufficient numbers of particles that statistical 
averaging is meaningful, but few enough so that there are seldom 
three or more bodies within the SOI at once. 

Consider a particular postion on the orbit of m2 
specified by its heliocentric distance, r2 • Encounters are 
possible with the fraction of m^ particles having periapses 
oriented such that the heliocentric distance to m^^ at the 
relative node of the orbits is nearly equal r2 . The 
heliocentric velocity components and relative velocity 
vector can be found from relations given by Opik ( 1976 , Interplane- 
tary Encounters ) . For a particular pair m^ and m^ , the patched 
conic model required transforming the orbit of m^^ from helio- 
centric to mass centered when m^ is at the sphere of influence of m2 • 
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The approach conditions of , relative position and velocity 
at the sphere of influence, define a hyperbolic encounter 
tractory of about m^. When is again at the sphere of 
influence on the outward leg of the trajectory, the orbit 
is "patched" back into a heliocentric conic orbit, which 
generally is different from the pre-encounter heliocentric 
orbit of The effect of the close encounter is dominantly 

to change the heliocentric velocity vector cf ra^ via ^ 
rotation of the relative velocity vector during the encounter. 
The angle through which the relative velocity vector 
rotates, v, depends on the impact parameter, b, relative speed, 
U, and the mass of the scattering body, m 2 , and is 


tan 


V 


m 


bU 


( 1 ) 


For any orientation of the orbit planes producing an 
intersection of the orbits, the relative speed is fixed. Hence, 
the post encounter heliocentric velocity is determined by 
the deflection angle v and the plane in which the rotation 
occurs, v/hich is the plane of the encounter hyperbolic 
trajectory relative to the scattering body, e^. The mean 
value of the post encounter velocity vector may be found by 
averaging over all deflection angles and relative inclinations 
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assuming a uniform areal distribution of impact parameters. 

The normalized probability distribution of deflection 

angles given the encounter speed U is given by Weidenschilling 


(1975, M 80, 145-153) . 




Where v 


■ / ' 0 YSvy 

W tr ‘win L 




1, 


), 


and v^. are deflection angles corresponding 
max min jr = 


to trajectories that just graze the surface of the scattering body 
and to trajectories that just enter the sphere of influence, res- 
pectively. The initial relative velocity vector expressed in a 
coordinate system with x radially outward from the Sun, y 
being orthogonal to x in the orbit plane of the scattering body 
and 2 being perpendicular to the scattering bodies' orbit 
plane is given by (see Figure 9) ; 


o 

U "tj 


[) ( CoSH-o po ) 

^ jj +• UKi rio £;nr. ^ 

T i_J S in ^ 




Geometry for rotation of relative velocity 
vector during gravitational scatter 



figure 9: 


■ ■ r 
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During an encourter in a fixed relative orbit plane, the 
velocity vector is rotated to a inal orientation given by 
= 0Q + V, hence the deflected relative velocity vector is 
found by replacing 3^ with 6^ in (3) . Then mean values are 


found from 










The final expressions are 


Ux = (F^Osr^ vir) + 4C tan (nv^ r tl) ) 

iT = {F^lva^,-XJ) + 4C tan F 2 Cm^r-tTri 

^ = ""zo '’1 ^2>'^ 

where Fj^(jTi 2 f>J) and F 2 depend only on the ratio 

^ , where is the two body escape speed from m 2 a tan 
^ U 

q _ 

j — -= ^ and C IS a parameter depending only on 

l/u + U 
y xo yo 


S^and is of order unity. 


I 
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To complote the averaging, we must average over all 
positions in the scattering bodies' orbit and over all 
orientafons of the relative node of the two orbits. Prior to 
carrying out these calculations, however, we shall numerically 
evaluate the initial averaged equations to determine the 
variations expected in further averaging. The averaged 
equations for the scattered velocity components were evaluated 
for several scattering body masses, initial orbits and orbit 
orientation and these results are summarized in Table 2. 

As can be seen from this table, the variation in mean 
scattered elements with position around the orbit is not large 
and for luost cases, the evaluation at quadrature provides 
a good approximation to the true mean when averaging over all 
positions in the orbit. Of major importance, however, is 
the variation with relative node, Afi. The encounter 
velocity is determined partially by the relative inclination 
which, for a given pair of orbital inclinations, is determined 
by Af2. The mean values of orbit velocity resulting from large 
numbers of randomly oriented individual encounters between 
bodies on given orbits is then found by averaging over all 
relative node orientations of the two bodies. 

These expressions for the mean scattered orbit e' ..!><■ :ts 
are being used to improve the velocity perturbation model in 
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the numorical simulation. During a timostop, tho number of 
encounters between the n^^ particles moving on orbit 1 and n^ 

particles moving on orbit 2 is calculated based on tho particle- 
in-a-box model along with tho new values of orbital elements 
of those particles that are scattered. Those scattered into 
new orbit element bins are removed from tho old bins and 
added to the now ones, while those that aren't deflected 
completely out of tho old bins are accounted for by 
recomputing tho moan value of orbit elements for the particles 
currently in the bin. By repeating this procedure for each 
pair of mass and orbit element bin, tho evolution of the 
orbit element distribution during the time step is found and 
a sequence of timestops gives tho evolution of the moan 
orbit element and hence tho moan impact speed distribution. 

The present version of the numerical model treats a single 
planetary zone so the semimajor axis is a constant and only 
variations in the eccentricity and inclination of planetesimals 
are modeled. Tho linking of different planetary zones, necessary 
to model the simultaneous evolution and growth of planets, 

(see Section II) may ultimately require a detailed treatment 
of the semimajor axis evolution. 
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TABLE 2: EFFECT OF POSITION IN ORBIT AND RELATIVE ORIENTATION 

ON MEAN DEFLECTION. 


Scattering Orbit 

Initial Scattered Orbit 

Final orbit when scattering 
position : 

Periapse 

Quadrature 

Apoapse 

Scattering Orbit 

Initial Scattered Orbit 

Final orbit when 

(A) Afi = 0 

Periapse 

Quadrature 

Apoapse 

(B) An = 90 
Periapse 
Quadrature 
Apoapse 

(C) A$1 = 180 
Periapse 
Quadrature 
Apoapse 


a(AU) 

e 

i (deg) 

.000 

.010 

o 

• 

o 

.000 

.050 

5.0 


.999 

.058 

4.94 

.998 

.059 

4.93 

.999 

.059 

4.94 

1.00 

.05 

4.9 

1.00 

.05 

5.0 


0.85 

.12 

5.2 

0.84 

.20 

5.3 

0.85 

.23 

5.3 

1.01 

.06 

11.8 

1.01 

.05 

11.8 

1.01 

.04 

11.8 

1.03 

.08 

14.5 

1.03 

.06 

14.5 

1.02 

.05 

14.5 
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VALIDATION OF MASS-SHIFTING ALGORITHM 


VALIDATION OF MASS-SHIFTING ALGORITHM 


12 

The large number of bodies (typically 10 ) involved in 

simulations of planet growth clearly makes it impossible to 
follow the mass evolution of individual objects. The usual 
approach in such problems is to treat the mass distribution 
as a continuum, and set up a generalized integro-dif f erential 
coagulation equation for its time variation (G.V. Pechernikova 
V.S. Safronov, and E.V. Zvyagina, Sov. Astron ., 20 , 346, 1976) 
However, the number ical solution of such equations requires 
very large amounts of computer time, even when the 
coagulation probability can be expressed as a simple function 
of the masses only. Our use of realistic mechanical proper- 
ties with a variety of collision outcomes, coupled with the 
simultaneous velocity evolution, makes this approach 
intractable. Instead, we use an approximate technique 
which conserves mass, but does not precisely account for 
numbers of particles. The basic method, described briefly 
by Greenberg et al. (1978) , will be reported here in greater 
detail, followed by a description of work in progress to test 
and improve its accuracy. 

(1) The mass distribution is divided into discrete size 
bins (usually spanning by a factor of two in mass) . Inter- 
actions between bins (collisions and gravitational stirring) 
are computed as if all particles in a bin have the mean mass 


of that bin. 
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(2) During a time stop, the number of bodies removed 
from a bin by accretion or catastrophic disruption in 
interactions with all other bins, and with itself, are 
calculated from particle-in-a-box statistics. The total 
net gain or loss of mass of the remaining bodies due to 
accretion or erosion is computed, and used to calculate the 
change in the mean mass in that bin. 

(3) Suppose that at the beginning of a ti..iostep, 
the total mass M in a bin is known. Wo have 

M = N in 

where N is the number of particles in the bin, and in 
is the mean mass. The change in mass during a timestep is 

dM = Ndiii + iiidN. 

The first term is evaluated in a straightforward manner, 
since the change in m has been computed, in the second term, 
the change in the number of bodies, dN, has tv/o components. 
One is due to those directly removed when accreted 
by other bodies or fragmented in catastrophic collisions. 

The other accounts for bodies which are assumed to lie near 
the boundaries of the bin, and can be "pushed" into another 
bin by a small change in mass. 

We assume that the number of particles in a bin is 
distributed uniformly in log M: 
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dN/d log m = C = constant across the bin. 


= / 


This gives ^ = J g— d log m = C log (m2/nij^) , 


where ~ lower and upper mass boundaries of the 

bin (Figure 1 ) . 

Then m = J’mdN{m) = C (m2 - nij^) and 

m = M/N = (m2 - mj^)/log (m2/m^) . 

Suppose, for example, that the mean change in mass per 
particle due to accretion and erosion in one time.step is 
Am. For simplicity, we assume that each particle in the 
bin had its mass multiplLed by the factor (1 + Am/m); this gives 
a uniform shift in log m across the bin. The mass 
is now uniformly distributed between m^^d + Am/m) and 
m3 = m^d + Am/in) (Figure 2 ). The mass pushed into the 
next bin is 


AM - J mdN (m) = C (m3 - m2) = Cm2Am/m. 

This mass is subtracted from the first bin and added 
to the next. The remaining mass in the first bin is no longer 
uniformly distributed, and must be adjusted before the 
next timestcp. After all mass changes have been computed, 
so the total mass in a bin is M', the number of particles 
is reset to N' = M'/ni/ and in the next timestep, we use the 
value C = N'/log . 


log m 




m. 


schematic of mass transfer algorithm. 


For 


clarity, only one bin is shown initially occupied. 


FIGURE 1 
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Since this algorithm is only approximate, its degree 
of accuracy is of some concern. An exact solution is not 
required, but we certainly want our results to give a 
qualitatively correct picture of planet growth within our 
initial assumptions, and hopefully to produce quantitative 
results within an order of magnitude of "exact" solutions. 

We are cautiously optimistic that this is indeed the case, 
for the following reasons: 

Analytical calculations using computed velocities and 
collision rates show that it is physically possible for the 
largest bodies to grow in the times elapsed; they do not 
appear to be artifacts of the program. 

We find that for the growth of the largest bodies, the 
first term in the expression for dM (Ndm) dominates. 

This term, due to accretion of smaller bodies, does not 
depend on the mass shifting algorithm. 

The program has been tested with a range of bin widths 
and timesteps, with little effect on the results. 

We have tried variations in the mass-transfer algorithm. 
For example, C has also been computed from the average number 
of particles in adjacent bins, to approximate a continuous 
distribution. All versions tested to date give qualitatively 
similar results. The time to produce a largest body of a 
given size typically varies by about a factor of two. 
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The program logic has been applied to some simple 
coagulation pi'obloms for which analytic and numerical 
solutions exist (Scott, 1968). In tested cases, it has 
computed mass distributions which are generally well within 
an order of magnitude of the exact solutions, provided the mean 
mass of the distribution has not evolved by more than an 
order of magnitude. This condition holds in all planet- 
growth simulations run to date. Additional tests are 
planned, and we are experimenting with other mass-transfer 
algorithms in hopes of further improvement. 

In extending our simulations to produce objects of 
planetary sise, we can avoid further dependence on such 
approximations. Once the largest bodies become sufficiently 
few in number, they can be detached from the size distribution 
and treated as individuals, rather than as part of a bin 
population. This "hybrid" approach will be incorporated into 
the program during the current year. 
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INCLINATION-TYPE RESONANCE 


IN A DISSIPATIVE MEDIUM 
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INCLIN7iTI0N-TYPE RESONANCE IN A DISSIPATIVE MEDIUM 

Earlier work on orbital resonances in dissipative media 
(Greenberg, 1978) led to exciting new results relevant to 
planet-building and to the structure of rings and of the 
asteroid belt. Past work focused on eccentricity-type 
resonances, the simplest to treat mathematically. Now we are 
looking at inclination resonances, which are generally more 
difficult to treat because the relevant equations are non-linear. 
We next describe current progress in treating this problem: The 

governing equations have been solved analytically. We plan in 
the immediate future to use this new solution as a basis for 
including a dissipative medium in the analysis. 

Consider perturbation of an interior satellite by a 
more massive perturber near the 2:1 resonance. Assume 
m << m' << M and i' = e = e' = 0. The resonance variable 
(critical argument of the disturbing function) is 

iji = 4A' - 2X - 20 . 



(The following analysis may not require the assumption that 
i' =e' =e=0, as long as conditions are such that there are 
no other important resonance variables nearby. See Greenberg, 
MNRAS 165, 1973.) 
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The disturbing function takes the following form after 
eliminating short-period terms and higher order terms in i: 

„ Gm'Pi^ cos '!i 1 Gm' .2 , ( + 1) 3 GM , , , ,2 .2 

R - — p ^ - 2 .. b3/2 - j -- JjCr/a) i 

(+1 ) 

where F and are both functions of a ^ a/a'. See 

3/ 2 

Greenberg (MNRAS 165, 1973) for details. Note F ^ 0.2 and 
^3/2 ^ Changing notation: 

R = F^Gm' (i^/a) cos J (i^/a) 

where F^(a) 0.13 

F2(a) % 1 


Also let Fj^' = 9F^/3a and F^' h 3F2/9a. 

Now 

• • • 

= 4n ' - 2n - ?e - 2fi 

where dots are time derivatives, c is mean longitude at epoch, 
and other elements are standard notation. Evaluate rhs using 
standard perturbation equations: 


’ = _ 2ria 9R na i 9R 

^ GM 3a GM 2 ^ 

2 


2na 

GM 


' Gm' (i^/a) cos i(/ - F^ Gm' (i^/a^) cos ip 
- F^' Gm' (i^/a) + F^Gm' (i^/a^) 


Q 9 9 4 

+ ~ GMJ.R i /a 
4 4 ' 


na 1 
GM 2 


1 2FTGm' (i 

2 1 a 


i 3 ^’^^^2^ i 
) cos * - 2 P 2 Gm' (-) - 2 — I— (|) 
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E = -2 (m'/M) ~ ^’jL^ 

+2 (m'/M) (^" 2 '^ “ ^2^ 

9 T /„2 . 2. .2 

-j J^ (R /a ) nx 

2 

+ (m'/M) Fj^ni cos ip 

- (m'/M) F^ni^ 

3 „2 / 2, .2 

- j (J^R /a )ni 

' _ na 3R 
^ GMi 9i 

= 2 (m'/M) F^ n cos ip - 2 (m'/M) F^n - — {J^^/'X^)n 


• • • . 2 

Assume ^ >> e because e cc x . 
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Also ov'aliuito n and i 


n 


n^a 

GM" \ 



na_ 

GMi 3i2 



To simplify define 
/ 


q 

4(m'/M) 

+ 3 (J 

2i-Vn 

di 

11 'dt 



n 

n/n' 



,t yi 

elds 



It 

di 

= 4 - (2 - 

q) n - 

4 (m ' 

da 

- - 2(m'/M) 


cos 1 || 

dn 

di 

= - 6(m'/M) 


cos i!' 


To si iify even further, assume that n as it appears in 
coefficients is nearly constant {n^) . That assumptio is 
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generally valid in this type of analysis, but should be 
checked carefully before final acceptance of results (cf. 
Greenberg ^ 78, 1973). Next let 

N “ 4 - (2 - q) n 
I =_ 1= 

C n - 2Fj^ n(m'/M) 


Now the governing equations take the form: 


^ = N + 2C cos ii< 
di 


dl 

dr 


= 2CI sin ij) 


= -3 (2 - q) Cln^ sin ijj 


Removing t dependence: 


Iy = - l,5(2-q)C 


N - Nq = -1.5(2 - q) Cn^Cl - 1^) 


dij; ^ Nq - 1.5(2 -q) Cn^d Iq) ^ ^ cos il> 

dl ^ ^ T- „ j „ .1. T 


2CI sin li' 


sin ij) 


Let (J) = -I cos ip 
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Then 

^ = I sin 'f' ^ ij; 

d(j, _ Nq - 1.5(2 - q) CnQ(I - 1^) 
dl “ 


which yields 


N + 1.5(2 - 

/f — ^ 

■ ‘3) Cn^lQ 3 „ x2 , 

2C 

0 ^ Q COilETZ- 


where the constant is evaluated by setting (f) = -1^ cos and 

T = T 

-0* 

We thus have a solution for the trajectory in (n, i, i|j) 
space. This is a closed trajectory. In order to include 
dissipative effects, such as imposed changes in n due to 
gas drag or tides or in i due to drag, we intend in the 
immediate future to apply a variation of parameters approach 
such as that used by Greenberg (M 7£, 1973) . 

Next consider the case in which i' ^ 0. In this case a 
resonance occurs where the resonance variable 

is slow varying. Let us assume that precession of and ' is 
such that this resonance is well separated from the one 
discussed earlier, i.e., that this ^ is the only slow-varying 
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argument. In this case, the equations of motion can be put 
in the form; 


^ = 4n' - 2n + (m'/M) n’ (i'/i) F cos i}; 

(m’/M) n'i'F sin \p 

6n (m'/M) ii'F sin \jj 

Since i' and F are nearly constant, these equations have 
the same form as the e-type resonance studied by Greenberg 
( Icarus 33 , 1978), except that e is replaced by i. Hence, 
an i damping medium in this case will have the same effect as 
the e damping medium had in that case, namely a secular 
variation in semi-major axis. The effect will be somewhat 
weaker due to the extra small factor (i') which appears here. 
We intend in the immediate future to consider the implications 
of these results in solar system evolution. These are some 
examples of our continuing technical work on resonance 
effects in the early solar system. 


di 

dt 


dn 

dt 
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THEORY OF ROTATIONS 

Work was initiated on understanding and modeling 
the collisional evolution of planetesimal rotation in con- 
junction with the population size frequency evolution. 

The focus of this year's research was rotational evolution in 
collisional fragmentation models, particularly with application 
to asteroids. Analytic models for equilibrium rotation states 
of asteroids have been developed by Harris ( Icarus in press 
1979) . The analytic rotation theory is applicable only 
to populations distributed according to a power law size- 
frequency distribution - a significant restriction in light 
of the plausible planetesimal, as well as observed asteroid, 
distributions. Large gravitationally bound asteroids are 
predicted to nave rotation rates independent of size, but they 
rotate considerably slower than is predicted using our best 
estimates of collisional parameters combined with the Harris 
theory: 3 hours vs. the 9 hour observed mean rotation 

period. Hence, for the asteroid belt, either (i) our 
estimates of parameters are incorrect, (ii) collisional 
relaxation has not yet been reached in the present belt, 
or (iii) the restriction to power law populations precludes 
model results from being applicable to asteroids. 

Option (i) is possible given that our knowledge of 
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impact physics at high speeds involving large bodies is 
quite limited. However, to reconcile theory and observation 
would require extreme values for some parameters, such as 
converting all the input collisional kinetic energy into 
ejecta kinetic energy. Optivon (ii) implies that there 
has been only a modest collisional evolution of asteroids, 
i.e., the asteroid belt was not substantially more massive 
initially than it is today. This result is consistent with 
results from collisional evolution modeling. Option (iii) 
is probably the most plausible, since the best determination 
of the asteroid size-frequency distribution indicates 
that the belt differs significantly from a power law distri- 
bution. Another shortcoming of the model is that the 
rotational models are applicable basically to a single 
generation population whereas most asteroids are multi- 
generation bodies, being fragments from collisions of larger 
bodies. The rotational aistribution of a fragmental popula- 
tion probably is determined in part by the spin rates of the 
fragments from the catastrophic disrUi’tion of large bodies, 
whereas in the theory developed so far, it is assumed that 
fragments rotate at the same rate on the parent body. But, 
this is unrealistic since the kinetic energy that goes into 
rotation is influenced more by the total kinetic energy of the 
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collisions than by the pre-impact rotation speed of the 
parent body. 

According to the model, for small bodies where material 
strength dominates over gravitational binding, the mean 
equilibrium rotation period is predicted to vary directly 
with the diameter of the asteroid - smaller bodies should 
rotate faster. Observations of asteroid rotations do 
suggest that small asteroids tend to rotate faster, although 
further observations are required to test the prediction of 
a linear dependence between rotation period and asteroid 
diameter. 

Further work is planned for the current year on 
extending models of collisional rotation. Primarily, it is 
planned to treat numerically the simultaneous evolution of 
the size-frequency distribution and the rotation distribution. 
By treating the problem numerically, not only can the 
evolution of both distributions be followed, but also the 
restrj.ction to power law distributions will no longer be re- 
quired and arbitrary population distributions can be studied. 
The rotational evolution model will be applicable to the 
rotational distribution for accreting planetesimals and 
should be applicable to studying the rotational properties of 
the growing planets. Improved rotational models may permit us 
to distinguish whether or not the observed rotation periods 
of large asteroids are primordial or the product of collisions. 
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ORIGIN OF ASTEROIDS; VELOCITY PUMPING BY LARGE JUPITER- 
SCATTERED PLANETESIMALS 

The present asteroid zone is unusual for its low 
mass density and high mean collision speed ('vS km/sec) j 
neither condition could have existed at the time asteroids 
formed. With the present mass density, intervals longer 
than the age of the solar system would be required to accrete 
Ceres-sized bodies, while a mean collision speed of 5 km/sec 
produces collisional destruction rather than accretion for 
asteroidal size bodies. However, the low mass density may 
be a consequence of the high collision speed, because 
collisional grinding will produce an asteroid-like mass 
distribution on a timescale less than the age of the solar 
system as shown by the collisional evolution models of 
Davis et aJ. (1979) , provided most of the mass is initially 
in bodies less than a few hundred km diameter. However, it 
has not been shown that velocity pumping proceeded mass 
depletion. For example, the same process that increased 
collision speeds may have also caused a mass loss, or conceiveably 
mass depletion could have occurred by a totally different process 
prior to velocity pumping. Jupiter is usually invoked as 
the agent £'v these effects. 

Var ous investigators, e.g., Safronov (1972, NASA TT F-677) 
and Weidenschilling (1974, Icarus 22, 426-435) have argued 
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that impacts of Jupiter-scattcred planetesima.ls (JSP) on 
early asteroids were responsible for the low mass density and 
high relative velocities. .Alternatively, Heppenheimer (1979, 

LSC-X Abs. 531-533) and Ward (1978, B.A.A.S . 2, 591) 
have suggested that Jovian secular resonances swept through 
the inner solar system as nebular gas dissipated and augmented 
planetesimal eccentricities and inclinations. Presumably 
the asteroid zone was affected to some extent by this latter 
hypothesized mechanism because there are still secular 
resonances in the belt today. However, let us consider first the 
the effects of JSP’s. Davis ^ (1979, see Appendix A. 4) 

argue that asteroids are too weak to have their velocities changed 
appreciably by collisions; they will be destroyed by large 
collisions before there are enough small impacts to signifi- 
cantly change their orbits. We suggested that gravitational 
stirring by large Jupiter-scattered planetesimals could have 
effectively pumped up asteroid velocities on a timescale com- 
parable with the elimination of JSP by Jupiter. In this 
scenario, a growing proto- Jupiter was accompanied by a size 
distribution of other bodies in tis zone, some of which were 
gravitationally deflected into the asteroid zone. 

To test the JSP scattering hypothesis, we have carried 
out numerical simulations using a version of our planet 
forming program that was modified to simulate asteroid collisional 
evolution as part of a separately funded program to study 
asteroids. The velocity pumping model is a 
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preliminary version of the algorithm described in Appendix B.4. 
The increase in mean asteroid collision speed with limt= is shown 
in Figure 1 fc- the case of an earth-sized body in a Jupiter- 
crossing orbit, but having its perihelion deep within the 
asteroidal zone. In this case, gravitational stirring oy the 
massive JSP increases the mear* asteroid collision speed from 
200 m/sec to 5 km/sec in x 10 years; an interval com- 
parable to the mean lifetime of bodies on Jupiter-crossing 
orbits. Subsequent collisional evolution among asteroids would 
reduce the initial asteroid population to the present observed 
pcpulati.on, provided that there were few additional original 
ast'^roid larger than 500 km diameter. 

The advantages of this JSP scenario for the early 
evolution of asteroidal velocities are: (1) there is no 

problem in eliminating the JSP, since the same agent (Jupiter) 
that brings them into the asteroid belt also removes them 
from the belt. (2) The stirring is accomplished by 
gravitational encounters rather than by collisions. (3) Large 
asteroids, such as Pallas with its 35° inclined orbit, are 
stirred as effectively as small bodies. (4) The reduction of 
mass in the asteroidal zone may be accomplished either by 
collisions with an accompanying population of smaller JSP 
or by collisions among the asteroids at their enhanced 
velosities, depending on the original mass distributions 
of both population components. (5) The existence of one or 


c 
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more Earth-sized JSP is a plausible accompaniment to the 
format 'on of the outer planets (although its necessity 
remains to bo demonstrated) and the timescale for their 
evolution is reasonable. 

Possible limitations of this scenario include: (1) 

There is no ready explanation of the Kirkwood gaps; they must 
be formed by a different process than the velocity pumping 
mechanism. (2) The degree of gravitational stirri-g depends 
on the random orbital evolution of the few largest Jupiter- 
scattered bodies, and therefore on the statistics o small 
numbers. More rigorous calculations are being planned to 
show that such an event, together with subsequent collisional 
evolution amoi\g the asteroids themselves, can produce the 
observed distributions of eccentricities and inclinations, 
as well as their mean values. 
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"Formative Conditions Controlling Structure of Planetary 
Highlands." An Abstract submitted to the Lunar Highlands 
Conference by William K. Hartmann 
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FORMATIVE CONDITIONS CONTROLLING STRUCTURE OF 
PLANETARY HIGHLANDS. William K. Hartmann, Planetary Science 
In.stituto, Tucson, Arizona 85719. 

Crater densities on highlands (or the most heavily cratered 
surfaces) of Moon, Mercury, Mars, Phobos, and Dcimos reach re- 
markably similar values (Fig. 1) , relative to the wide scatter 
of densities in other j^rovinces. This coincidence may bo tolling 
us something about highland structure, and suggests some mecha- 
nism that defined an envelope .of crater densities rarely cxcoodod 
on planetary landscapes. This process is likely to be a satu- 
ration or steady state mechanism, implying repeated pulveriza- 
tion of highlands, producing deep mogaregoliths . 

Reconstruction of early accrctionary cratering (Fig. 2) 
can bo done invoicing (1) present-day observed flux rates; (2) 
observed decay of meteorite flux from 4 Gy ago to present, based 
on dated lunar samples; and (3) flux rate necessary to build 
Earth and Moon in time intervals compatible with meteorite for- 
mation intervals. The impact flux in the "missing" years from 
4.5 to 4.0 Gy ago ranged about 10^ to 10' times the present flux 
rate, had an average 20 to 45 my half-life, caused enormous rates 
of megaregolith production, and docs not appear to re quire any 
anomalous "spike" of catastrophic cratering around 4.0 Gy ago. 

Relative timing between magma ocean cooling, and decay rate' 
of cratering was very important in controlling present day ob- 
served character of planetai-y highlands, since it controlled 
outcome of competition between preservation of intact igneous 
rock and conversion of igneous roclcs into massive breccias. 
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FIGURE 1. Crater densities on "highlands" provinces on several 
planets lie near values about 32 x those found among large cra- 
ters in lunar maria. Clustering is tight relative to diversity 
of values in other provinces . 
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accretion of Earth in 0.1 to 150 my 
and 

accretion of Moon in 0.03 to 100 my 
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FIGURE 2. Schematic reconstruction of accretionary cratering 
without invoking a "spike" of catastrophic cratering 4 Gy ago. 
Present rate is observed (lov/er right) and projected back to 4 
Gy ago using known impact rates determined from lunar dating. 
Additional datum is impact rate required to form Earth and Moon 
during their estimated formation intervals. Connection in 
"missing interval" from 4.5 to 4.0 Gy gives plausible planetesi 
mal sweepup half-lives of the order 20 to 45 my. 
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"Early Runaway Growth of Planets." An Abstract 
submitted to the Division of Planetary Science 
meeting of the American Astronomical Society 
in St. Louis, 1979, by Richard Greenberg 
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Early Run.'iwav Crov.’lh of rianot !j. RICHARD 
GREENBERG, Plain.- tary Sci. liiat. (Tuc.ion) Exten.sion of 
our numerical simulation of planet cz-owth beyond the 
stage previously rej^orted (Greenberg cj_ Icarus . 35 , 

1, 1978) indicates that the largest body in an accretion 
zone grows so muc)i faster than smaller bodies tnat, by 
the time it is a few thousand km ,in diameter (only 
-5x10^ yr after collisional accretion begins) , it is 
detached from the continuous size distribution of 
planetcsimals. As it becom.es detached, but not before, 
it reaches a size such that a significant fraction of 
the zone's mass resides in tiie largest body. This 
scenario contradicts assum.ptions in analytical m.odels 
of planet growth (Safronov, NASA TT F-677, 1972) that 
the largest bodies contain most of the m.ass through- 
out most of the growth period and that they all rem.ain 
part of the continuous portion of the size distribution. 
We are searching for material and system parameters 
that might yield evolution consistent with Safronov's 
assumption. Our model does give rel-ative velocities 
of planetesimals consistent with Safronov's analytical 
results ; The velocities are on the order of the 
escape velocities of those bodies which dominate 
the continuous portion of the size distribution. We 
have sim.ulated the growth of bodies nearly 8000 km 
in diameter. 


